Drift instability and tunneling of lattice solitons 
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We derive an analytic formula for the lateral dynamics of solitons in a general inhomogeneous 
nonlinear media, and show that it can be valid over tens of diffraction lengths. In particular, we 
show that solitons centered at a lattice maximum can be "mathematically unstable" but "physically 
stable" . We also derive an analytic upper bound for the critical velocity for tunneling, which is valid 
even when the standard Peierls-Nabarro potential approach fails. 
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PACS numbers: 42.65 Jx, 42.65 Tg, 03.75 Lm 



Solitons have been thoroughly studied in diverse fields 
of physics such as nonlinear optics, BEC, plasma and wa- 
ter waves. By now, the stability and dynamics of solitons 
in homogeneous media are well understood. The possibil- 
ity to manufacture transparent materials with spatially 
varying, high contrast dielectric properties raises new 
questions regarding stability and dynamics of solitons in 
inhomogeneous media. In particular, while in a homoge- 
neous medium the solitons can freely move sideways, in 
an inhomogeneous medium the loss of translation invari- 
ance affects the lateral movement of solitons. The prob- 
lem of lateral movement of lattice solitons is interesting 
theoretically and important for applications such as all 
optical switching and quantum information science. It 
was studied analytically, numerically and experimentally 
for various media and lattices, see e.g., 0, 0, 0, 0, H, [fif- 
However, each of these studies considered a specific non- 
linearity, lattice type and dimension. 

In this Letter, we provide, apparently for the first time, 
a unified theory for the mobility of lattice solitons which 
is valid for any nonlinearity, lattice type and dimension. 
We show that soliton mobility is intrinsically related to 
soliton stability, two key properties that so far were stud- 
ied separately. This relation enables us to compute ana- 
lytically the rate of drift of solitons initially centered near 
a lattice maximum, and the restoring force that the lat- 
tice exerts upon solitons initially centered near a lattice 
minimum. In the latter case, our approach provides an 
upper bound for the critical velocity for tunneling, which 
is valid even when the standard Peierls-Nabarro potential 
approach cannot be applied. 

All solitons centered at a lattice maximum are "math- 
ematically unstable" as they drift towards the nearest 
lattice minimum 0,0, HQ- However, the ability to com- 
pute the magnitude of the drift rate allows us to identify 
cases in which the drift rate is so small so that the soliton 
is "physically stable", i.e., the drift instability does not 
develop over the propagation distance of the experiment. 
This observation explains why in some experiments soli- 
tons centered at a lattice maximum were observed to be 
stable Q. 

Consider the dimensionless nonlinear Schrodinger 



equation (NLS) with a linear lattice 

iA z (z, x) = -V 2 A - F (\A\ 2 ) A + V(Nx)A, (1) 

where x = (xi, . . . , Xd) and V is a linear lattice/potential 
with a characteristic length scale or period 1/N. This 
equation describes propagation in a media with Kerr 
nonlinearity as well as quintic, cubic-quintic, satu- 
rated/photorefractive nonlinearities etc. The variable z 
denotes the propagation coordinate in nonlinear optics 
and time coordinate in BEC. Eq. ([1]) has soliton solu- 
tions A = e IAtz u(x), where u is the solution of 



V 2 w(x) + F{\u\ 2 )u - y(iVx)w — flu = 0. 



(2) 



It is well known that a necessary condition for stability 
is the slope (VK) condition ^ > where V = f u 2 dx. is 
the soliton power. Violation of the slope condition leads 
to a width instability, i.e., small perturbations can lead 
to large changes of the soliton width, which in some cases 
result in collapse 0, IE> @] • 

If the soliton is centered at a lattice maximum, it 
can become unstable even if the slope condition is sat- 
isfied 0, 0, S, S 0|- Indeed, there is a second con- 
dition for stability, the spectral condition, which states 
that for A — ue l ^ z to be stable, the number of neg- 
ative eigenvalues of the operator L 



(V) 



V^(x) — F(u 2 ) — 2u 2 F' should be at most one more 
than the number of negative eigenvalues of the operator 



(V) _ 



-V 2 + M + ^(x) - F(u 2 ) 



L 

A simpler version of the spectral condition was derived 
in [4j for the case u > 0. Recall that in a homogeneous 
medium, has d zero eigenvalues Ag^ — = with 



corresponding eigenfunctions f!j V °' = J^-, where Q = 
u (v=o) ig the so i u tion of (2]) with V = 0. The potential 
V breaks the translation symmetry of the medium. As a 
result, XqO can split into d different values. The spectral 

condition is violated if and only if at least one XqO attains 
a negative value [4|. 

The spectral condition can be derived from the 
following linear stability analysis Q. Let A = 
e vz (u(x) + h(z, x)) where h(z,x) is a small perturba- 
tion. Since the instability due to violation of the spectral 
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condition originates only from the eigenfunctions /,■ (x) 



of iS^l which correspond to negative A y • , we can rewrite 
the perturbation h as 

h(z, x) = Cl e°^(/f } + igf } ) + c_ie- 
where c±i are constants and 
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Since L^/f > 



Aq -■ and is positive definite, 



-CyAg 



(V) 



(V) f (V)\ 



f 00 ,00* 



>0. (3) 



Therefore, fij is real (i.e., instability) when Aq^ < and 
imaginary (i.e., stability) when Xqj > 0. 

The effect of a lattice on A ^ was studied in e.g., [4, 
HL S S > where it was shown that if the soliton is cen- 
tered at a lattice minimum (maximum), then Xqj be- 
comes positive (negative), hence the spectral condition is 
satisfied^ violated) . 

In WMM 

it was observed numerically that violation 
of the spectral condition results in a drift instability, i.e., 
the center of mass (COM) of the beam in the Xj coordi- 
nate, defined as xj(z) = J Xj\A\ 2 /V, drifts away from its 
initial location xj(0) near the lattice maximum. So far, 
however, the relation between the spectral condition and 
the drift instability has not been established analytically. 

odd i, 



To do that, we note that since and are 

1 2 

Xj(z) = —(xj,\u(x;fj,) + h(z,x)\ ) 



(4) 



= B ( Cl e n > z +c_ie 



-n,- 



where B — 2(xj,ufj)/V is constant. Thus, 



(5) 



x 3 {z) = VjXjiz). 

Relation ^ shows that a failure to satisfy the spectral 
condition (Xqj < ) leads to a drift instability. More- 
over, the magnitude of flj = |CVAqj |= determines the 
drift rate away from the lattice maximum in the Xj di- 
rection. 

A simpler expression for Cy can be obtained for a weak 
lattice (V <C /i) and a power nonlinearity F = \A\p~ 1 



In this case, /, 

-(V=0) 



(V) 
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(V=0) 



- % and I,™ * 
From the Pohozaev identities it follows that (-P^-, §0-) = 

\ ox j ' axj * 

jg(Q,Q) where S — 1,2 ^p^ +d ■ In addition, if we mul- 
tiply L- ^w — J^L by X jQ and integrate in parts we get 
( L -^§§-> Wj) =\(Q,Q)- Substituting in © gives 



-(V=0) 



n 2 



M,00 



(6) 



Hence, for a weak lattice, the dependence of the drift 
rate VLj on the lattice period N is only through its ef- 
fect on Aq . The approximation (|6|) can be generalized 
for different nonlinearities and to lattices which are not 
weak. For example, in the case of narrow lattice, Q, 2 is 
given by ((6|) with fj, replaced by a + V(xo) where Xo is 
the location of the soliton peak Q . 

We solve Eq. ([1]) numerically for F = \A\ 2 , d — 1 and 



V{x) = V cos(2ttNx), 



(7) 



with the initial condition A(0, x) — u{x — S), where u(x) 
is the solution of (J2|) centered at x = 0. Therefore, x{z = 
0) = 5 (since d = 1, we can suppress the index j). For a 
small shift S <C 1, we can rewrite A(Q, x) = u(x) + h(0, x) 
where h(0,x) = -6% + 0(S 2 ). Since h(0,x) = {a + 
c-i)/ (y) +?;(ci-c_i) 5 ( y ), then ci-c_i = and by g]), 



x(z) — Scoshttz. 



(8) 



In our simulations we observe that that the COM 
evolves according to ((SJ, see Fig. Qla), and therefore, 
calculate the drift rate numerically by finding the best 
fitting f2. We fix Vq = 0.1 and /i = 4.5 and vary N. As ex- 
pected for a soliton centered at a lattice maximum 0, [f| , 



Aq^ < for all values of N, and Xq" ' vanishes in the lim- 
its N — > (narrow solitons) and N — * oo (wide solitons), 
see Fig. [TJb) . In Fig. QJc) we confirm that the numer- 
ically computed drift rate £1 is in excellent agreement 
with Eq. §3§ and also with the approximation ©. Ac- 
cordingly, each value of A^ is attained at two different 
values of N for which the drift rates are nearly identi- 
cal. Indeed, in Fig. [TJa) we see that the drift rate of 
the COM when N = 0.2 and N = 1.437, both for which 
Aq = —0.0454, is the same over more than 3 orders of 
magnitude and 40 diffraction lengths. In Fig. [IJd) we 
repeat these simulations with a stronger lattice (Vo = 2). 
In this case, the numerically computed drift rate is in 
excellent agreement with the one predicted by Eq. ([3]). 
The approximation © is very accurate only for narrow 
(N <C 1) and wide (N ^> 1) solitons. Indeed, although 
the lattice oscillations are not small, for narrow and wide 
solitons, the effect of a mean-zero lattice is weak, hence 
the deviation of f^ from t^L is small 0, @. Although 

for solitons of N = 0(1) width the deviation of /j ^ from 
^ is not small, the approximation ([5]) is, at most, 10% 
inaccurate. 

In Fig. [He) we fix N = 1 and Vo =0.1 and vary /i. 
As in Fig. [2(b), since the soliton is centered at a lattice 
maximum, Aq^ < for all values of /i, and X^ van- 
ishes in the two limits /j, — ► (wide solitons) and fi — > oo 
(narrow solitons). In Fig.QJf) we see that fi is monotoni- 
cally increasing in /i, and that the numerically calculated 
drift rate is in excellent agreement with the analytical 
prediction ([3]) and also with its approximation 
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We emphasize that despite the similarity of the depen- is nowadays frequently used in optics and BEC [l(| • 

In [11] it was proved that PIM converges only if l} V%> 



dence of Aq ' onJV and fi (see Fig. QJb) and Fig. [He)), 
the dependence of O on TV and /i is completely differ- 
ent in the narrow-beam limit. Indeed, for narrow beams 



A, 



(V) 



-AN 
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2d 2 V | 



jV 2 d 2 V I 
u dx 2 I 



„ =0 @ so that by ©, ft 2 w Q 2 , 



dx 2 \x=0' 



Hence, il vanishes for a fixed fi and 
TV — > (Fig. HJc)) but approaches n narroul = 2.8 for 
a fixed TV and \i -> oo [Fig. [JJf)]. 

In order to show that our results are also valid in higher 
dimensions, we solve Eq. flj in a d — 2 setting with 



V(x, y) — (cos 2 {2-kx) + cos 2 (27ry)) 



with Vo = 5, and find the numerical drift rate to be in 
excellent agreement with Eq. ([3]), see Fig. [Tig). Remark- 
ably, although the lattice is strong, the numerical drift 
rate is also in excellent agreement with the approxima- 
tion ([6]) in which is shifted by Vn/2, the mean of V. 




FIG. 1: (Color online) (a) The dynamics of the COM for 
Vo = 0.1, S = 1CT 4 and the lattice © for TV = 0.2 (dashed 



red line), TV = 1.437 (dotted blue line), and the 
prediction ([8j) with Q ~ 0.52 (black solid line). The 
indistinguishable, (b) Ag as a function of TV. (c) 
Q as a function of TV. The analytical prediction ([3]) 
line) and its approximation © (dashed red line) 
indistinguishable, (d) Same as (c) for Vo = 2. (e) 
function of fi. (f) Same as (c), but as a function 
Same as (f) for d = 2 and the lattice ([9]) with Vo = 



analytical 
3 lines are 
Drift rate 
(solid blue 
are nearly 

(V) 


of /J,. 
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The analytical relation ([5]), together with © or ©, 
enable us to estimate the distance at which a soliton ini- 
tially centered near a lattice maximum will deviate sig- 
nificantly from its initial location. In particular, if the 
initial shift S and/or drift rate f2 are sufficiently small, 
then this "mathematically unstable" soliton can remain 
"near" its initial location over the propagation distance of 
the experiment, i.e., be "physically stable". This obser- 
vation can explain the experimental results of [3J, where 
solitons centered at a lattice maximum did not drift over 
ss 18 diffraction lengths. 

The relation between the sign and magnitude of the 
perturbed near-zero eigenvalues {Aq^}^ =1 and the drift 
instability appears to be universal. Indeed, we now show 
that it also occurs in numerical calculation of soliton pro- 
files using Petviashvili's iterations method (PIM), which 



has at most one negative eigenvalue, which is a spectral 
condition similar to the one for the stability of NLS soli- 
tons. Accordingly, PIM is not expected to converge for 
solitons centered at lattice maxima. 

We solve Eq. © with F = | u\ 2 using PIM with the 
initial guess — u(x — 5), where u(x) is the solu- 
tion of centered at a maximum of the lattice (O with 
Vo = 0.1. Similarly to the dynamics of NLS solitons cen- 
tered slightly off a lattice maximum, the COM of 

(data not shown), 



(9) evolves according to x(m) 



where u( m ) is the solution in the mth iteration. Thus, 
we conclude that when the spectral condition for PIM is 
violated, the method does not converge because the iter- 
ative solution drifts away from the lattice maximum. In 
that sense, the analogy between the dynamics (in z) of 
NLS solitons and of vS™"' (in m) is further demonstrated, 
since in both cases, violation of the spectral condition 
leads to a drift instability. We also compute the expo- 
nential drift rate numerically for various combinations 

of TV and fi and observe that £1 2 = — Dy ^Aq^' (TV, fi) / /ij 

where the constant Dy depends on Vq but is independent 
of fi and TV, see Fig. [UJa). Interestingly, the scaling of 
in Aq V '' > and in /i is different from ©, yet in both cases 
fl depends on TV only through Aq^'. 

Although in [11] it was proved that the iterations 
should diverge for solitons centered at a lattice max- 
imum, in several studies these iterations did "con- 
verge" [1, 0, [12] ■ To explain this apparent inconsistency, 



(m) 



as a function of 



in Fig. [2Kb) we plot max x \u 
for TV = 0.2, [i = 2 and the lattice with V = 0.1 
and = e~ x , and observe that the iterations con- 
verge (i.e., max x \u^ m ^ — u\ < 10~ 13 ) after « 40 iterations. 
However, if we continue the iterations, a significant drift 
of the COM occurs around m f=a 2000. To understand 
this "post-convergence" drift, we note that in this exam- 
ple, S -0.085 and n S 0.018. Since the seed of 
the drift is roundoff error, then x(m = 0) = O(10 -16 ). 
Indeed, io- 16 e 018 - 2000 =0(1). This example of a nu- 
merical iterative solution which theoretically should di- 
verge yet in practice converges is thus analogous to the 
mathematically unstable yet physically stable NLS soli- 
tons discussed earlier. 
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FIG. 2: (Color online) (a) Drift rate fi as a function of X Q V ^ /fj, 
for the solution of © using PIM. (b) Maximal error (solid 
line) and COM (dashed red line) in the mth iteration. 

We now consider solitons centered near a lattice min- 



4 



imum. Since A^ > 0, the spectral condition is satis- 
fied. Hence, these solitons are stable under small lateral 
perturbations. Indeed, relation © shows that small lat- 
eral perturbations would lead to oscillations around the 
lattice maximum, while relation ([3]) shows that the mag- 
nitude of determines the strength of the restoring 
force. For example, consider a soliton centered at a lat- 
tice minimum which is launched at an angle 9j between 
the Xj and z axes. Such an angle corresponds to an ini- 



tial transverse velocity of uqj = = 0) 
©-(gD, the COM evolves according to 



Xj{z) = u ,jsin(|Oj|2)/|Oj 



tan#j. By 



(10) 



Thus, as Ag j , 



hence IfL-l, increase, the maximal devi- 



ation of the COM from the lattice minimum becomes 
smaller, implying stronger lateral stability. 

Eq. (TIT))) gives an accurate description of the dynamics 
for small velocities. However, for non-small velocities, as 
the soliton propagates sideways, the attraction towards 
the lattice minimum decreases, an effect which is not cap- 
tured by Eq. (fTO)) . To see that, we solve Eq. © with 
d = 2, F — \A\ 2 and the lattice © with V = 0.5. The 
initial condition is A(0,x,y) = u(x,y)e^ v ° x+Voy)/2 , i.e., 



a soliton centered at a lattice minimum x„ 



(0,0) 



with initial velocity in the direction of the nearest lat- 
tice maximum at x ma2; = (0.25,0.25). Indeed, for small 
initial velocities, the agreement between the dynamics 
and Eq. (fl"0|) is excellent, see Fig. [3(a). For higher ve- 
locities, the COM initially evolves according to Eq. (fT0|) 
but deviates from it as it approaches the lattice maxi- 
mum, see Fig. [3(b). For a sufficiently large initial ve- 
locity, the soliton can "tunnel" beyond the nearest lat- 
tice maximum. The critical velocity for tunneling Vg r 
is the one for which the transverse velocity x(z) van- 
ishes at the lattice maximum. The upper limit |v§ r | < 



r: 



=i i 2 ( x « 



Xmm)? can be derived from 



th y Zjj=l l"JI v«TfM*x '~mm/j 

Eq. ITU)) . In the case of Fig. [31(b), this bound gives |vg r | < 
1.22, an over-estimate of w 45% over |vg r | w 0.485V2 . 



(a) 
max 




FIG. 3: (Color online) Dynamics of xi(= X2) (solid) and 
theoretical prediction (|10[) (dashes) in a bulk medium with 
F = \A\ 2 , the lattice © with Vo = 0.5, fi = 35 and (a) 
v = (0.2,0.2), (b) v = (0.485,0.485). (c) Power of solitons 
centered at a lattice maximum (solid) and minimum (dashes) . 

The standard formula for vjf , based on the Peierls- 
Nabarro potential (PNP) approach is | vg 1 " | = 

y/AATl/V where AH is the difference in the Hamiltoni- 
ans of equal-power solitons centered at a lattice minimum 




FIG. 4: (Color online) Same as Fig. for F = \A\ 2 -0.02|yl| 4 , 
fi = 2.5, Vo = 1 and (a) v = (0.0125,0.0125), (b) v = 
(0.01905,0.01905). 

and maximum, respectively. In order to apply the PNP 
approach, the power of the soliton centered at a lattice 
maximum should be equal to that of a soliton centered at 
a lattice minimum. For a two-dimensional Kerr medium 
(F = \A\ 2 ), however, such "soliton pairs" do not exist, 
since the power of all solitons centered at a lattice max- 
imum is below that of all solitons centered at a lattice 
minimum, see Fig. [3(c). Therefore, one cannot use the 
PNP approach, and the upper bound provides the 
only analytic estimate of |vq' |. 

Finally, we solve Eq. © for a cubic-quintic nonlinear- 
ity and the lattice ©. As in the Kerr case, for small 
initial velocities, the agreement between the numerics 
and Eq. (|10p is excellent over many diffraction lengths 
[Fig. Ufa)], while for higher velocities the COM deviates 
from Eq. (JTDJ) as the soliton approaches the lattice maxi- 
mum [Fig. E(b)] . In this case the PNP approach is appli- 
cable [see Fig. [3(c)] and yields |vg r | ~ 0.027. The value 
of the critical velocity obtained numerically is within 1% 
of the PNP prediction [see Fig. [4(b)]. To the best of our 
knowledge, this is the first demonstration of quantitative 
agreement of the PNP approach with numerical results 
for d — 2. The research of G.F. and Y.S. was partially 
supported by BSF grant no. 2006-262. 
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